
% This routine is used to calculate the strain spectra from MMP profiles at
% all the three MMP locations (1, 3, and 4). 
%
% Author: Ritabrata Thakur, 23 November 2021
%
% Last used: 30 March 2022

clear

MMP_no_array = {1,3,4};

for jj=1:length(MMP_no_array)  
% This automatically chooses the MMP number from the above array
MMP_number = MMP_no_array{jj}; 

clearvars -except jj MMP_no_array MMP_number 
fprintf('Running MMP_number %d\n',MMP_number);

%% Loading data
% this is depth and time range where the velocity is good (not NaN)
% for MMP1
if (MMP_number == 1)
 MMP_load=load('iwap06_MP1_WHOI112.mat');
 depth_start = 43;
 depth_end   = 700;
 lat = 25.5;   %North
 lon = 195;  %East
 start_time = 4;
 end_time = 980;
 inertial_start = 25.05; % 90% of 27.84 hrs
 inertial_end   = 30.6;  % 110% of 27.84 hrs
end
% for MMP3
if (MMP_number == 3)
 MMP_load=load('iwap06_MP3_WHOI114.mat');
 depth_start = 44;
 depth_end   = 702;
 lat = 28.9;   %North
 lon = 196.5;  %East
 start_time = 4;
 end_time = 980;
 inertial_start = 22.34; % 90% of 24.83 hrs
 inertial_end   = 27.31; % 110% of 24.83 hrs
end
% for MMP4
if (MMP_number == 4)
 MMP_load=load('iwap06_MP4_101.mat');
 depth_start = 44;
 depth_end   = 700;
 lat = 30.1;   %North
 lon = 197.1;  %East
 start_time = 4;
 end_time = 538; %this is different from the others because this is what is available data
 inertial_start = 21.52; % 90% of 23.92 hrs
 inertial_end   = 26.31; % 110% of 23.92 hrs
end

% reading the raw variables. For strain we do not need velocities
mmp.temp = MMP_load.mp.t;
mmp.salt = MMP_load.mp.s;
mmp.depth= MMP_load.mp.z; 
% the following array of one hour is from Savage et al dataset.
mmp.time_array  = squeeze(MMP_load.mp.datenum1h);

% time range to calculate the spectra from. Not used currently as the
% original time array has some issue. Using the mmp.time_array which was
% probably made up for use in Savage et al

%% Brunt Vaisala
mmp.Z = -mmp.depth(depth_start:depth_end);
mmp.p = gsw_p_from_z(squeeze(mmp.Z),lat);
mmp.depth_array_short = -mmp.Z; 

[SA, ~] = gsw_SA_from_SP(mmp.salt(depth_start:depth_end,:),mmp.p,lon,lat);  %converts practical salinity
%to absolute salinity which is used by the N2 routine. The flag of "1" if
%the measurement is in ocean. If it is used in land, set it to 0. do not
%put anything. it is automatically set
CT = gsw_CT_from_t(SA,mmp.temp(depth_start:depth_end,:),mmp.p); %converts in-situ temperature to
%conservative temperature as is required in GSW toolbox

clear N2 NN
% Brunt Vaisala frequency
for i=1:length(CT(1,:))
 [N2,p_mid] = gsw_Nsquared(SA(:,i),CT(:,i),mmp.p,lat);
 NN(i,:)=N2;  %this is the total Brunt-vaisala frequency squared stored
 % at a single location. It is probably best to continue with one location
 % as averaging over different locations would mean we would have to
 % calculate the wkb factor for each depth differently as the
 % stratification would vary according to location
end

% removing the values where there is unstable stratification
for i=1:length(NN(:,1))
    for j=1:length(NN(1,:))
      if(NN(i,j) <= 0)
      NN(i,j) = NaN;
      end
    end
end

%sqrt of the N^2
BV = sqrt(NN);

%% Strain calculation

%interpolating to original depth levels. As the original units are in
%pressure (p), so interpolating it to that level would amount to
%interpolating to original depth. The original depths are equally-spaced
%with a separation of 2m. Hence we can do an usual FFT

for i=1:length(NN(:,1))
 BV_interp(i,:) = interp1(p_mid, BV(i,:), mmp.p, 'linear', 'extrap');
end

for i=1:length(NN(:,1))
 NN_interp(i,:) = interp1(p_mid, NN(i,:), mmp.p, 'linear', 'extrap');
end

BV_interp_time_mean   = mean(BV_interp, 'omitnan');
BV_interp_time_mean   = naninterp(BV_interp_time_mean);
NN_interp_time_mean   = mean(NN_interp, 'omitnan');
NN_interp_time_mean   = naninterp(NN_interp_time_mean);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Strain calculation
for i=1:length(mmp.depth_array_short)
  strain_mmp(:,i) = (NN_interp(:,i) - NN_interp_time_mean(i))./NN_interp_time_mean(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Interpolate over NaNs (this is basically the naninterp function)
strain_mmp(isnan(strain_mmp)) = interp1(find(~isnan(strain_mmp)), strain_mmp(~isnan(strain_mmp)), find(isnan(strain_mmp)), 'spline', 'extrap');

%% Filtering the strain data into various frequency bands. I am having to
% not use some times because the values are NaNs. This was initially not
% done but filtering the strain into different frequency bands on 30 March.
% We need not filter BV and can directly filter strain and calculate the spectra.

start_time = 2;
end_time   = length(strain_mmp(:,1))-1;
for i=1:length(strain_mmp(1,:))
  %strain_mmp_3_to_5hr(:,i)  = bandpass(strain_mmp(4:end-1,i),[1/5.20 1/2.80], 1); 
  %strain_mmp_6_to_8hr(:,i)  = bandpass(strain_mmp(4:end-1,i),[1/8.20 1/5.80], 1); 
  strain_mmp_highpass(:,i)  = highpass(strain_mmp(start_time:end_time,i),1/11.5, 1); 
  strain_mmp_lowpass(:,i)   = lowpass(strain_mmp(start_time:end_time,i),1/11.5, 1); 
  strain_mmp_sem_diur(:,i)  = bandpass(strain_mmp(start_time:end_time,i),[1/13.5 1/11.5], 1); 
  strain_mmp_inertial(:,i)  = bandpass(strain_mmp(start_time:end_time,i),[1/inertial_end 1/inertial_start], 1); 
end
strain_mmp = strain_mmp(start_time:end_time,:);

%% WKB scaling.
N_0_mean = 0.0052; %used in literature = 3cph
mmp.z_stretched(1) = mmp.depth_array_short(1);
for ii=2:length(mmp.depth_array_short)
 mmp.z_stretched(ii) = mmp.z_stretched(1) + trapz(mmp.depth_array_short(1:ii), squeeze(BV_interp_time_mean(1:ii)))./N_0_mean; 
end

% this array is needed to interpolate the velocities from unequally-spaced
% stretched coordinates to equally-spaced
mmp.z_stretched_equal_spaced = linspace(mmp.z_stretched(1), mmp.z_stretched(end), length(mmp.z_stretched(1:end)));
N = length(mmp.z_stretched_equal_spaced);
win = hann(N);

% Number of FFT points for the Fourier transform
NFFT = N; 

%Fs=2*pi/(mmp.z_stretched_equal_spaced(2)-mmp.z_stretched_equal_spaced(1)); %removed 2pi on 04Nov21
%following early literature (Gregg et al 1991: "varieties of fully resolved spectra of vertical shear")
% this 2pi is also removed from shear and strain
delta_z = 1/(mmp.z_stretched_equal_spaced(2)-mmp.z_stretched_equal_spaced(1)); 
mmp.waveno = 0:delta_z/NFFT:delta_z/2; % we can do this because the depth
% array is equally spaced albeit in stretched coordinates

% getting the time-averaged psd from the following. I am also getting
% "psdx_array" which is dependent on time too. This has been output to
% understand the time evolution of the psd

%% Strain PSD and saving
clear wave_psd_avg_strain_mmp wave_psd_avg_mmp_strain_inertial wave_psd_avg_mmp_strain_6_to_8 wave_psd_avg_mmp_strain_3_to_5
% equally-spaced depth array is used because I have already interpolated
% the Brunt Vaisala to the original depth array "mmp.p" array in which 
% the depths are separated by 2m.
%%%%%%%%%%%%%%%%
% [psdx_array_strain_mmp_3_to_5,mmp.wave_psd_avg_mmp_strain_3_to_5]     = waveno_only_spectra_strain_ybp(strain_mmp_3_to_5hr,mmp.z_stretched_equal_spaced,NFFT);
% [psdx_array_strain_mmp_6_to_8,mmp.wave_psd_avg_mmp_strain_6_to_8]     = waveno_only_spectra_strain_ybp(strain_mmp_6_to_8hr',mmp.z_stretched_equal_spaced,NFFT);
% [psdx_array_strain_mmp_inertial,mmp.wave_psd_avg_mmp_strain_inertial] = waveno_only_spectra_strain_ybp(strain_mmp_inertial',mmp.z_stretched_equal_spaced,NFFT);
[~, mmp.wave_psd_avg_strain]                       = waveno_only_spectra_strain_ybp(strain_mmp,mmp.z_stretched_equal_spaced,NFFT);
[~, mmp.wave_psd_avg_strain_inertial]              = waveno_only_spectra_strain_ybp(strain_mmp_inertial,mmp.z_stretched_equal_spaced,NFFT);
[~, mmp.wave_psd_avg_strain_highpass]              = waveno_only_spectra_strain_ybp(strain_mmp_highpass,mmp.z_stretched_equal_spaced,NFFT);
[~, mmp.wave_psd_avg_strain_lowpass]               = waveno_only_spectra_strain_ybp(strain_mmp_lowpass,mmp.z_stretched_equal_spaced,NFFT);
[~, mmp.wave_psd_avg_strain_semdiur]               = waveno_only_spectra_strain_ybp(strain_mmp_sem_diur,mmp.z_stretched_equal_spaced,NFFT);
%%%%%%%%%%%%%%%%


mmp.wave_psd_avg_strain              = mmp.wave_psd_avg_strain./(NFFT*delta_z);
mmp.wave_psd_avg_strain_inertial     = mmp.wave_psd_avg_strain_inertial./(NFFT*delta_z);
mmp.wave_psd_avg_strain_highpass     = mmp.wave_psd_avg_strain_highpass./(NFFT*delta_z);
mmp.wave_psd_avg_strain_lowpass      = mmp.wave_psd_avg_strain_lowpass./(NFFT*delta_z);
mmp.wave_psd_avg_strain_semdiur      = mmp.wave_psd_avg_strain_semdiur./(NFFT*delta_z);

mmp_psd_strain.waveno     = mmp.waveno;
mmp_psd_strain.tot_strain = mmp.wave_psd_avg_strain;
mmp_psd_strain.inertial   = mmp.wave_psd_avg_strain_inertial;
mmp_psd_strain.highpass   = mmp.wave_psd_avg_strain_highpass;
mmp_psd_strain.lowpass    = mmp.wave_psd_avg_strain_lowpass;
mmp_psd_strain.semdiur    = mmp.wave_psd_avg_strain_semdiur;
mmp_psd_strain.comment = 'The MMP depth range is from 83-1348m for the whole period from 25 April 2006';

basedir_strain = '/home/';
savename = [basedir_strain 'mmp' sprintf('%d', MMP_number) '/psd_strain.mat'];
if exist(savename, 'file')
    delete(savename);
end
save(savename, 'mmp_psd_strain');

end




